www.gusucode.com > Hilbert-Huang Transformmatlab程序,及word版程序详单,这种算法用于机械行业故障诊断 > Hilbert-Huang Transformmatlab程序,及word版程序详单,这种算法用于机械行业故障诊断/源程序/源程序/hspec.m
% HSPEC: Hilbert Amplitude Spectrum % % [S,freq]=hspec(imf,N); % % S - Time-frequency-amplitude matrix % Columns are indexed in time, rows in frequency, values are amplitudes % % freq- instantaneous frequencies of each component % % imf - Matrix of intrinsic mode functions (each as a row) % % N - Number of frequency cells % % See: Huang et al, Royal Society Proceedings on Math, Physical, % and Engineering Sciences, vol. 454, no. 1971, pp. 903-995, % 8 March 1998 % % Remark: the graphical representation is the Hilbert Energy Spectrum % that is: 20 * log ( S * S ) % % Author: Ivan Magrin-Chagnolleau <ivan@ieee.org> % function [S,freq] = hspec(imf,N); L = size(imf,1); % Number of components in the decomposition %------------------------------------------------------------------------- % loop for on each component S = []; % Matrix which will contain the time-frequency-amplitude representation clear x z m p freq x = imf'; % now each column is a component z = hilbert(x); % analytic signal m = abs(z); % module of z p = angle(z); % phase of z for i = 1:L freq(:,i) = instfreq(z(:,i)); % instantaneous frequency % if the function instfreq is not available... % p(:,i) = unwrap(p(:,i)); % unwrap phase % freq(:,i) = abs(diff(p(:,i))); % derivative of the phase and absolute value % to have always positive frequencies ceilfreq(:,i) = ceil(freq(:,i)*N); % to have integer values - also do a smoothing given the number % of frequency cells for j = 1:length(x)-2 S(ceilfreq(j,i),j+1) = m(j+1,i); end end eps = 0.00001; % to avoid zero values before the log E = 20 * log ( S.^2 + eps ); % Hilbert energy spectrum % plot S figure; t=1:length(x); % time sample f=t/length(x)*0.5; % normalized frequency imagesc(t,f,E); % !!! I am not sure it is the best way to visualize it !!! colorbar; set(gca,'YDir','normal'); xlabel('Time Sample'); ylabel('Normalized Frequency'); return